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Abstract 

Tight-frame, a generalization of orthogonal wavelets, has been used successfully in various problems in image processing, 
including inpainting, impulse noise removal, super-resolution image restoration, etc. Segmentation is the process of identifying 
object outlines within images. There are quite a few efficient algorithms for segmentation that depend on the variational approach 
and the partial differential equation (PDE) modeling. In this paper, we propose to apply the tight-frame approach to automatically 
identify tube-like structures such as blood vessels in Magnetic Resonance Angiography (MRA) images. Our method iteratively 
refines a region that encloses the possible boundary or surface of the vessels. In each iteration, we apply the tight-frame algorithm 
to denoise and smooth the possible boundary and sharpen the region. We prove the convergence of our algorithm. Numerical 
experiments on real 2D/3D MRA images demonstrate that our method is very efficient with convergence usually within a few 
iterations, and it outperforms existing PDE and variational methods as it can extract more tubular objects and fine details in the 
images. 

Index Terms 

Tight-frame, thresholding, image segmentation, wavelet transform. 

I. Introduction 

Segmentation problem of branching tubular objects in 2D and 3D images arises in many applications, for examples, extracting 
roads in aerial photography, and anatomical surfaces of blood vessels in medical images. In this paper, we are concerned with 
identifying tube-like structures in Magnetic Resonance Angiography (MRA) images. Because of the necessity to obtain as much 
fine details as possible in real time, automatic, robust and efficient methods are needed. However, due to the low contrast, 
intensity inhomogeneity, partial occlusions and intersections, and the presence of noise and possible blur in the images, this 
segmentation problem is very challenging. 

There are many different approaches for vessel segmentation, see for example lfl6l . Il20l . Il23l . 11261 . PP . 041 . 11421 and the 
extended reviews BJ, (30 1. In the following, we concentrate on two approaches that are related to our approach: the active 
contour approach and the partial differential equation (PDE) approach. For the vessel segmentation algorithms that are based 
on deformable models [35], because the explicit deformable model representation is usually impractical, level set techniques to 
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evolve a deformable model have been introduced, and they provide implicit representation of the deformable model. However, 
the level set segmentation approach is computationally more expensive as it needs to cover the entire domain of interest, 
which is generally one dimension higher than the original one. Interested readers can refer to recent literature on the level set 
segmentation strategy for tubular structures |26], ||29l , BOl . 

Based on the curve evolution techniques, Mumford-Shah functional and level sets, a new model for active contours to detect 
objects in a given image was proposed in IT6l . Unlike the classical active contour models, this model does not depend on the 
gradient of the image. In [38], a generalization of the active contour without edges model was proposed for object detection 
using logic operations. However, both the active contour model [16] and the logic framework [38 1 are not suitable for detecting 



tubular structures with low contrast and intensity inhomogeneity, see numerical results in Section IV 

In ||23l , a geometric deformable model for the segmentation of tubular-like structures was proposed. The model is charac- 
terized mainly by two components by using a suitable diffusion tensor: the mean curvature flow and the directionality of the 
tubular structures. The major advantage of this technique is the ability to segment twisted, convoluted and occluded structures 
without user interactions; and it can follow the branching of different layers, from thinner to larger structures. The dependence 
on the grid resolution chosen to solve the discretized PDE model is still an open problem. The authors in [23 1 have also 
applied a variant of the proposed PDE model to the challenging problem of composed segmentation in [24|. 

Besides the methods above, there were some initial work by using wavelets or tight-frames to do texture classification and 
segmentation Oj, J43j. The tight-frame approach is a versatile and effective tool for many different applications in image 
processing, see [4], (7), fi"2l . 1131 . [39|, [41], 11441 . There are many kinds of tight-frame systems, such as those from framelets 
ll36l . contourlets |[T9l and curvelets iflOl . ATI , etc. Recently, the authors in [20] proposed to combine the tight-frame image 
restoration model of [8[ and the total variation based segmentation model of Q, |fl5l , |fl6l to do segmentation. Their approach 
results in a minimization problem. In this paper, we also derive a segmentation algorithm that uses the tight-frames. However 
our method is not based on minimizing any variational model and hence it is different from the method in 11201 . In fact, our 
algorithm iteratively updates the set of possible boundary pixels to change the given image into a binary image. Like the 
method in 1231 . our method also has the ability to segment twisted, convoluted and occluded structures. We will show that our 
algorithm is convergent. In fact, numerical results show that it converges within 10 iterations for 2D as well as 3D Magnetic 
Resonance Angiography (MRA) images. We will see that our method can extract much more details from the given image 
than the method in lfl6l . [20 1 and [23 1. We remark that a preliminary version of our segmentation algorithm has been given 
in our proceeding paper (9). The main contributions in this paper are: (1) a simpler strategy to initialize and refine the regions 
enclosing the possible boundary of the vessels, (2) the new strategy leads to a simple proof of convergence and easier choice 
of parameters, (3) a different tight-frame with better directionally selective property is used to obtain more details in the image, 
(4) a new 3D result with higher noise and more complicated vessel is added. 

The rest of the paper is organized as follows. In Section [nl we recall some basic facts about tight-frame and tight-frame 



algorithms. Our segmentation algorithm is given in Section III In Section IV we test our algorithm on various real 2D and 



3D MRA images and compare it with those representative algorithms from different approaches: PDE-based [23], tight-frame 
ll20l and active contour models |H51 . Conclusions are given in Section [V] 
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II. Tight-Frame Algorithm 



In this section, we briefly introduce the tight-frame algorithm which are based on tight-frame transforms. All tight-frame 
transforms A have a very important property, the "perfect reconstruction property": A T A = I, the identity transform, see 
11361 . Unlike the wavelets, in general, AA T 7^ T, For theories of framelets and tight-frame transforms, we refer the readers 
to ifTTI for more details. In order to apply the tight-frame algorithm, one only needs to know the filters corresponding to the 
framelets in the tight-frame. In the followings, we give two examples of tight-frames: the piecewise linear B-spline tight-frame 
11181 and the dual-tree complex wavelet tight-frame |37|. 

The filters in the piecewise linear B-spline tight-frame are: 



ho = *[1, 2, 1], ht = ^[1, 0, -1], h 2 = i[-l, 2, -1], 



(1) 



see 11361 . The tight-frame coefficients of any given vector v corresponding to filter hi can be obtained by convolving hi with 
v. In matrix terms, we can construct, for each filter, its corresponding filter matrix which is just the Toeplitz matrix with 
diagonals given by the filter coefficients, e.g. Hq — |tridiag[l, 2, 1]. Then the ID tight-frame forward transform is given by 

H 



A = 



Hi 
H 2 



(2) 



To apply the tight-frame transform onto v is equivalent to computing Av, and i7;v gives the tight-frame coefficients corre- 
sponding to the filter hi, i = 1,2,3. 

The d-dimensional piecewise linear B-spline tight-frame is constructed by tensor products from the ID tight-frame above, 
see |2T1 . For example, in 2D, there are nine filters given by = hf <£> hj for i,j = 1,2,3, where hi is given in (fTJ. For 
any 2D image f, the tight-frame coefficients with respect to fty are obtained by convolving hij with /. The corresponding 
forward transform A will be a stack of nine block-Toeplitz-Toeplitz-block matrices (cf. Q). The tight-frame coefficients are 
given by the matrix- vector product Af, where f = vec(/) denotes the vector obtained by concatenating the columns of /. 

Dual-tree complex wavelet transform (DCWT) was firstly introduced by Kingsbury 11321 . [33). Apart from having the usual 
perfect reconstruction property, shift-invariance property and linear complexity, it also has the nice directionally selective 
property at ±15°, ±45°, ±75°. The idea is to use two different sets of filters: one gives the real part of the transform and the 
other gives the imaginary part. Let {go, gi} and {ho, hi} denote the two different sets of orthonormal filters, where go, ho are 
the low pass filters and g%, hi are the high pass filters. Let the square matrix A 3i h, denote the 2D separable wavelet transform 
implemented using g. L along the rows and hj along the columns, and define Ah i f l .,A gig .,Af lig . similarly, where i,j = 1,2. 
Then the forward transform for the 2-dimensional DCWT is represented by 
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We see that the 2-dimensional DCWT requires four different wavelet transforms in parallel. We refer the readers to fl4l . ifTTl . 
Il28l . 11371 and the references therein for more details and for the implementation of 3-dimensional DCWT. In practice, the 
DCWT are implemented by using different sets of filters for the first level and the remaining levels, see [371 ■ The filters and 
the Matlab code for DCWT can be obtained from lfl4"l . 

The tight-frame algorithms, as given in J4), Q, 0, Q, iTHl . are of the following generic form: 

f('+l) = W(fW), (3) 

f (i+i) = A T Tx{Af (t+ ^) 1 i = 1,2, (4) 

Here fW is an approximate solution at the i-th iteration, U is a problem-dependent operator, and 7a( ) is the soft-thresholding 
operator defined as follows. Given vectors v = [i>i,--- ,t>„] T andA = [Ai,--- ,A„] t ,7a(v) = [t\ t (fi), • • • ,t\ n (v n )] T , where 

J sga(v k )(\vk\ - Xk), if\v k \>Xk, 
tx k {Vk) = < (5) 
[ 0, if |« fc |<Afc. 

For how to choose Xk, see ll22l . 

We remark that Q performs a tight-frame denoising and smoothing on the image while Q performs a data-fitting according 
to the specific problem at hand. For high-resolution image reconstruction problem 1T31 . astronomical infra-red imaging J5], 
impulse noise removal 0, and inpainting problem |4|, the tight-frame algorithm has been shown to be convergent to 

a functional with the regularization term being the 1-norm of the tight-frame coefficients, see (7). 

III. Tight-Frame Based Algorithm for Segmentation 

The technology of Magnetic Resonance Angiography (MRA) imaging is based on detection of signals from flowing blood 
and suppression of signals from other static tissues, so that the blood vessels appear as high intensity regions in the image, 
see Fig. [TJa). The structures to be segmented are vessels of variable diameters which are close to each other. In general in 
medical images, speckle noise and weak edges make it difficult to identify the structures in the image. Fortunately, the MRA 
images contain some properties that can be exploited to derive a good segmentation algorithm. From Fig. [TJa), we see that 
the pixels near the boundary of the vessels are not exactly of one value, but they are in some range, whereas the values of the 
pixels in other parts are far from this range. Thus the main idea of our algorithm is to approximate this range accurately. We 
will obtain the range iteratively by a tight-frame algorithm. The main steps are as follows. Suppose in the beginning of the ith 
iteration, we are given an approximate image and a set A^) that contains all possible boundary pixels. Then we (i) use 
AW to estimate an appropriate range [aj,/3j] that contains the pixel values of possible boundary pixels; (ii) use the range to 
separate the image into three parts — those below the range (background pixels), inside the range (possible boundary pixels), 
and above the range (pixels in the vessels); (iii) denoise and smooth the inside part by the tight-frame algorithm to get a new 
image f^ +1 \ We stop when the image becomes binary and this happens within 10 iterations for the real MRA images we 
tested, see Table [1] in Section [IV] In the followings, we elaborate each of the steps. Without loss of generality, we assume that 
the given image / has dynamic range in [0, 1]. 

Initialization. To start the algorithm at i = 0, we need to define /l°>, the initial guess, and A(°), the initial set of possible 
boundary pixels. Naturally, we set /W — /, the given image. For A", since we do not have any knowledge of where the 
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boundary pixels are at the beginning, we identify them by using the gradient of /, i.e. we locate them as pixels where the 
gradient is bigger than a threshold e. More precisely, let £1 be the index set of all the pixels in the image, then we define 

A (0) = {j e Q | ||[V/]j||i > e}. (6) 

Here [V/]j is the discrete gradient of / at the jth pixel. Once /(°) and A^ ^ are defined, we can start the iteration. 

Step (i): computing the range [oti, Pi]. Given A^, we first compute the mean pixel value on A^: 

u (l) = V f (l) (7) 

» i A (i)i 2^ h ' V> 

1 1 ,-eAW 

where | • | denotes the cardinality of the set and /j- is the pixel value of pixel j in image /W. Then we compute the mean 
pixel values of the two sets separated by ^M; 



and 



1 

Ki.AW : 7->^}| {ieA(< ^^ i)} 

While /i^ reflects the mean energy of the set of possible boundary pixels, ix_ and fjX} reflect the mean energies of the pixels 
on the boundary closer to the background and closer to the vessels respectively. We define 

fM (i) +^ ) l l 

ai=max< ,0>, ftsmm — (10) 

Step (ii): thresholding the image into three parts. Using the range [oj,&] C [0,1], we can separate the image /W into 
three parts — those below, inside, and above the range, see Fig. |TJb) for i — 0. Since our aim is to create a binary image, we 
threshold those pixel values that are smaller than a,; to 0, those larger than to 1, and those in between, we stretch them 
between and 1 using a simple linear contrast stretch, see l27l . If there are no pixels in between on and then the threshold 
image is binary and the algorithm stops. More precisely, let 

Mi = max{/, (l) | a t < ff < ft, j G A«}, m, = min{/f 5 | a, < jf < € A«)}, (11) 



then we define 



0, if ff < at, 

fer' ^</f<ft, foraUjea (12) 

1, ifft</f°, 

Fig. [TJc) shows the threshold and stretched image from Fig. |TJb), where the yellow pixels are pixels we have classified not 

(j+i) (i+h) 

on the boundary, i.e. f - 2 = (signifying pixel j is in the background) or f - 2 — 1 (signifying pixel j is inside the 
vessel). The remaining pixels are remained to be classified and we denote the set by: 

A (i+1) ={j | 0</j i+ ' ) <l,jeO}. (13) 

Note that those pixels with values m, and M, are mapped to and 1 respectively and hence will not be in A( l+1 ). Next we 
denoise and smooth f( l+ h) on A^ +1 ^. 
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(a) 



(b) 



(c) 



Fig. 1. (a) Given MRA image, (b) Three parts of the given image (green-below, red-in between, and yellow-above), (c) Threshold and stretched image by 
(T5| (yellow pixels are with value or 1). 

Step (iii): tight-frame iteration. To denoise and smooth the image on A^ +1 ', we apply the tight-frame iteration (Wl) 



on A( l+1 ). More precisely, if j £ A^ +1 \ then we set = f^ +2 ^', otherwise, we use (4i to get To write it 

clearly, let f^+i) = vec(/^ + s- ) ), and be the diagonal matrix where the diagonal entry is 1 if the corresponding index 

is in hy +1 >, and otherwise. Then 



out 



f(*+i) = rj_ p(*+i))f(i+s) + p( l+ ^A T T\{Af {i+ ^ ] ). 



(14) 



By reordering the entries of the vector f( t+1 ) into columns, we obtain the image Note that the effect of (14i is to 

denoise and smooth the image on see 0. Since the pixel values of all pixels outside A(' i+1 ) are either or 1, the cost 



of the tight- frame transform in ( 14 1, such as the computation of ^4f ( l + 2 ) ? can be reduced significantly by taking advantage of 



the fact that all computations can be done only on pixels around A^ +1 '. 

Stopping criterion. We stop the iteration when all the pixels of /f'+s) are either of value or 1, or equivalently when 
AW = 0. For the binary image f^ + ^\ all the pixels with value are considered as background pixels and pixels with value 
1 constitute the tubular structures. 
Below we give the full algorithm and show that it always converges to a binary image. 

Algorithm 1: Tight-frame algorithm for segmentation 



1. Input: given image /. 

2. Set /(°) = / and A(°) by \6 1 

3. Do i = 0, 1, . . . , until stop 



(a) Compute by (10). 

(b) Compute f( l +^ by Q. 

(c) Stop if f( l+ h) is a binary image. 

(d) Compute A^ +1 ) by 



13 1. 



(e) Update f( i+ ^ to by 
Output: binary image 



14) 



Theorem 1. Our tight-frame algorithm will converge to a binary image. 



Proof: From (13 I, it suffices to prove that |A^| = at some finite step j > 0. By (111, if /(*+a) is not yet a binary 
image, then there will be at least one j <G such that /j 1 ' = Mi. By jTij , /j i+2 ' will be set to 1 and hence by ( joj , 
j £ A (i+1) . Hence |A (i+1) | < |A (i) |. Since |A (0) | is finite, there must exist some i such that |A W | =0. ■ 

We emphasize that the algorithm actually converges within 10 iterations for the 2D and 3D real images we have tested. 

Finally, let us estimate the computation cost of our method for a given image with n pixels. Since the costs of computing 
H®, l£\ and [a;,/3i] are all of 0(n) operations, see |t|-([io|i; and the cost of a tight-frame transform is also linear with 
respect to n (see e.g. (|2| where Hi are all tri-diagonal matrices), we see that the complexity of our algorithm is 0(n) per 
iteration. In fact, one can speed up the computation tremendously as all computations can be done only on pixels around AW 
and there is no need to carry out the computations in ft — though in the numerical tests, we did not optimize the code and we 
just carried out the computations in O. According to Table [i] after just 3 iterations, the set AW contains only a hundred pixels 
in the 2D case and just two thousand pixels in the 3D case. 

IV. Numerical Examples 

In this section, we try our tight-frame segmentation algorithm (Algorithm 1) on 2D and 3D real MRA images tested in 
ll23l and 11241 . Our algorithm is written in Matlab. The thresholding parameters used in (|5]l are all chosen to be = 0.1; 
and we choose e = 0.003 and 0.06 for 2D and respectively 3D images in |6]). To display the boundary and the surface in the 
binary images visually, we use the Matlab commands "contour" and "isosurface" for 2D and 3D images respectively. 
The tight-frame we used is the DCWT downloaded from 1 14 1 where all parameters in the code are chosen to be the default 
values. In particular, the number of wavelet levels used in the code is 4. The 2D images are tested in a MacBook with 2.4 GHz 
processor and 4GB RAM, while the 3D images, because of their sizes, are tested on a node with 120GB RAM in a PC-cluster. 
We compare our method with those representative algorithms from different approaches: PDE-based 1231 , tight-frame ||20l and 
active contour models fl6l . where the programs are provided by the authors. 

Example 1. The test image is a 182 x 182 MRA image of a carotid vascular system, see Fig. |2|a). The blood vessels contain 
regions with high and low intensities, including some very thin vessels in the middle with intensities as low as the intensity of 
the background. Intersections of partial structures even increase the difficulty of the segmentation. The segmentation results by 
different methods are given in Fig. [2] where the given image is overlaid so that we can compare the accuracy of the methods. 
Clearly, the results of Figs.[2|b) and (c) are not satisfactory since the vessels obtained in Fig. |2jb) are disconnected and some 
vessels in Fig. |2jc) are not detected. By comparing the parts inside the rectangles in Fig. |2|d) with those in Fig. |2je), we see 
that our method can extract clearer boundary than the method in 1231 . especially for handling the artifacts near the boundary. 
We also see that the noise in the image has been removed because of the denoising property of formula Q as explained in 
Section [n] Our method converges in 5 iterations and requires 0.64 seconds. The first column of Table [i] gives |AW| at each 
iteration. 



Example 2. The test image is a 256 x 256 MRA image of a kidney vascular system as shown in Fig. [3ja). This example 
shows the ability of our method to reconstruct structures which present small occlusions along the coherence direction. Our 
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(a) (b) (c) (d) (e) 

Fig. 2. Carotid vascular system segmentation, (a) Given image, (b), (c) and (d) Results by the methods in fl6|, |20| and |23| respectively, (e) Result of our 
method. 

TABLE I 

Cardinality of AW at each iteration of the four examples 



|A«| 


Example 1 
|n| = 11284 


Example 2 
|fi| = 66049 


Example 3 
|fi| = 8120601 


Example 4 
[fi| = 6000000 


i = 


1721 


9444 


137330 


152898 


i = 1 


354 


1943 


32760 


32064 


i = 2 


82 


464 


8795 


8565 


i = 3 


26 


133 


2475 


2391 


i = 4 


4 


32 


689 


650 


i = 5 





8 


189 


187 


i = 6 







56 


59 


i = 7 






15 


12 


i = 8 






3 


2 


i = 9 













method converges in 6 iterations with 0.78 seconds; and the second column in Table [ijgives |AW| at each iteration. The results 
by the method in lfl6l (Fig. [3jb)) and by the method in |20| (Fig. |3jc)) are not good since they can not recover the small 
occlusions along the coherence direction, while this can be done by our method and the method in (23], see Figs. |3jd) and 
(e). Furthermore, our method is better than the method in [23| by comparing the rectangular parts of Fig. |3jd) with those in 
Fig. [3|e), since our method can detect smoother edges; see Figs. [3jf)-(k) which are the results of zooming in the rectangular 
parts of Figs. |3jd) and (e) respectively. The results also show that our method is very effective in removing artifacts. 

Example 3. This is a 3D example where we extracted a volumetric data set of size 201 x 201 x 201 from a 436 x 436 x 540 CTA 
(Computed Tomographic Angiography) image of the kidney vasculature system, see Fig. |4fa). Because of different curvatures, 
diameters, bifurcations, and weak surfaces impaired by noise, it is hard to detect the tips of the thin vessels. Figs. |4jb) and 
(c) give the results by using the method in [24J and our proposed method respectively. The figure shows that our method can 
give much more details. A visual comparison of the given image with our result shows that almost all the vessels are correctly 
segmented. We note that our method converges in 9 iterations, see the third column in Table [I] 



(a) (b) (c) 




(f) (g) (h) 




(i) (j) (k) 

Fig. 3. Kidney vascular system segmentation, (a) Given image, (b), (c) and (d) Results by the methods in |16|, |20| and [23 1 respectively, (e) Results by 
our method, (f)-(k) are the zoomed-in parts of (d) and (e). 

Example 4. Our second 3D example is an MRA data set of a brain aneurysms (vessel wall dilatations). The 120 x 250 x 200 
volumetric data set has been extracted from a 120 x 448 x 540 MRA image of the brain-neck vasculature system, see Fig. |5ja). 
As in Example 3, the different curvatures, diameters and bifurcations of the vessels make it a difficult problem. In addition, the 
high noise makes the thin vessels hard to see even by the naked eyes. For this topological complex tubular structures, Figs. |5jb) 
and (c) give the results by the method in l23l and our method respectively. Obviously, our method can segment many more 
thin vessels, especially those corrupted by the high noise. The results reflect that our new method is very effective in handling 
the noise spread on the surface of the vessels. For this complicated example, our method converges in 9 iterations only, see 
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(a) 





(b) (c) 

Fig. 4. Segmentation of the kidney volume data set. (a) Given CTA image, (b) Result by the method in 1241 . (c) Result by our method. 

the fourth column of Table [I] Our segmentation technique can be used to compute vessel radii and other clinically useful 
measurements in case of aneurysms. In Fig. |5jc), one may argue that there are some small isolated points in the image and 
the surface of the vessels is not smooth. This can easily be remedied by smoothing our final binary image by the tight-frame 
formula Q one time before we show the image, since Q has the denoising property as we explained in Section [II] See Fig. 
|5Jd) for the denoised-and-smoothed image after one iteration of Q. 

V. Conclusions and Future Work 

In this paper, we introduced a new and efficient segmentation method based on the tight-frame approach. The numerical 
results demonstrate the ability of our method in segmenting tubular structures. The method can be implemented fast and give 
very accurate, smooth boundaries or surfaces. In addition, since the pixel values of more and more pixels will be set to either 




Fig. 5. Segmentation of the brain volume data set. (a) Given MRA image, (b) Result by the method in |23|. (c) Result by our method, (d) Result of our 
method after smoothing by |4j once. 
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or 1 during the iterations, by taking advantage of this, one can construct a sparse data structure to accelerate the method. 
Moreover, one can use different tight- frame systems such as those from contourlets, curvelets or steerable-wavelet fl9l . ifTTI . 
Il25l to capture more directions along the boundary. Though we have proved that our algorithm will always converge to a binary 
image, it will be interesting to see what functional the binary image is minimizing. The framework for proving convergence 
for tight-frame algorithms, as developed in |7|, may be useful here. These are the directions we will explore in the future. 
Acknowledgments.: This work was supported by HKRGC Grant CUHK 400510 and CUHK Direct Allocation Grant 2060408. 
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